
##load data
x.vars = c('non.white','hispanic','college.grad','male')
money.dat = read.csv('money_dat.csv')
attribution.dat = read.csv('attribution_dat.csv')
social.dat = read.csv('social_dat.csv')
physical.dat = read.csv('physical_dat.csv')
all.dat = read.csv('minimal_groups_dat.csv')


cat('\n \n Starting Experiment 2 Estimation \n \n') 

##outmatrix
outmat = matrix(nrow = 4, ncol = 11)
rownames(outmat) = c('allocation','attribution','social.perceptions','physical.perceptions')
colnames(outmat) = c('int.mean','int.sd','seg.mean','seg.sd','ate','se','ci.low','ci.high','p.value','cohens.d','N')

###allocation analysis
##RI analysis
money.Xs = as.matrix(money.dat[,x.vars])  ##get necessary x variables
##generate RI permutations and probabilities
money.ri.perms = genperms(money.dat$treatment.assignment,
                        blockvar = money.dat$Block.ID)
money.ri.probs = genprob(money.ri.perms)

##estimate treatment effect and inference
ate = estate(Y = money.dat$allocated.diff,
             Z= money.dat$treatment.assignment,
             prob = money.ri.probs,
             X = money.Xs)
money.Ys = genouts(Y = money.dat$allocated.diff,
                   Z= money.dat$treatment.assignment,
                   ate = 0)
##generate outcomes
distout <- gendist(money.Ys,money.ri.perms, prob=money.ri.probs, X = money.Xs)
money.Ys.ci = genouts(Y = money.dat$allocated.diff,
                      Z= money.dat$treatment.assignment,
                      ate = ate)
distout.ci <- gendist(money.Ys.ci,money.ri.perms, prob=money.ri.probs, X = money.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

##cohen's D for stadndardized effects
cohen = cohen.d(d = money.dat$allocated.diff,
                f = as.factor(money.dat$treatment.assignment))

##add outcomes to matrix
outmat['allocation','int.mean'] = mean(money.dat$allocated.diff[money.dat$treatment.assignment == 0])
outmat['allocation','int.sd'] = sd(money.dat$allocated.diff[money.dat$treatment.assignment == 0])
outmat['allocation','seg.mean'] = mean(money.dat$allocated.diff[money.dat$treatment.assignment == 1])
outmat['allocation','seg.sd'] = sd(money.dat$allocated.diff[money.dat$treatment.assignment == 1])
outmat['allocation','ate'] = ate
outmat['allocation','se'] = sd(distout.ci)
outmat['allocation','ci.low'] = quantile(distout.ci,.025)
outmat['allocation','ci.high'] = quantile(distout.ci,.975)
outmat['allocation','p.value'] = disp$greater.p.value
outmat['allocation','cohens.d'] = cohen$estimate
outmat['allocation','N'] = nrow(money.dat)


#################################
##now do same thing for all other outcomes
#################################

#####attribution analysis
attribution.Xs = as.matrix(attribution.dat[,x.vars])
attribution.ri.perms = genperms(attribution.dat$treatment.assignment,
                                blockvar = attribution.dat$Block.ID)
attribution.ri.probs = genprob(attribution.ri.perms)
ate = estate(Y = attribution.dat$attribution.diff,
             Z= attribution.dat$treatment.assignment,
             prob = attribution.ri.probs,
             X = attribution.Xs)
attribution.Ys = genouts(Y = attribution.dat$attribution.diff,
                         Z= attribution.dat$treatment.assignment,
                         ate = 0)
distout <- gendist(attribution.Ys,attribution.ri.perms, prob=attribution.ri.probs, X = attribution.Xs)
attribution.Ys.ci = genouts(Y = attribution.dat$attribution.diff,
                         Z= attribution.dat$treatment.assignment,
                         ate = ate)
distout.ci <- gendist(attribution.Ys.ci,attribution.ri.perms, prob=attribution.ri.probs, X = attribution.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = attribution.dat$attribution.diff,
                f = as.factor(attribution.dat$treatment.assignment))

outmat['attribution','int.mean'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 0])
outmat['attribution','int.sd'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 0])
outmat['attribution','seg.mean'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 1])
outmat['attribution','seg.sd'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 1])
outmat['attribution','ate'] = ate
outmat['attribution','se'] = sd(distout.ci)
outmat['attribution','ci.low'] = quantile(distout.ci,.025)
outmat['attribution','ci.high'] = quantile(distout.ci,.975)
outmat['attribution','p.value'] = disp$greater.p.value
outmat['attribution','cohens.d'] = cohen$estimate
outmat['attribution','N'] = nrow(attribution.dat)


####social analysis
diff.Xs = as.matrix(social.dat[,x.vars])
diff.ri.perms = genperms(social.dat$treatment.assignment,
                         blockvar = social.dat$Block.ID)
diff.ri.probs = genprob(diff.ri.perms)
ate = estate(Y = social.dat$social.diff.scale,
             Z= social.dat$treatment.assignment,
             prob = diff.ri.probs,
             X = diff.Xs)
diff.Ys = genouts(Y = social.dat$social.diff.scale,
                  Z= social.dat$treatment.assignment,
                  ate = 0)
distout <- gendist(diff.Ys,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
diff.Ys.ci = genouts(Y = social.dat$social.diff.scale,
                  Z= social.dat$treatment.assignment,
                  ate = ate)
distout.ci <- gendist(diff.Ys.ci,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = social.dat$social.diff.scale,
                f = as.factor(social.dat$treatment.assignment))

outmat['social.perceptions','int.mean'] = mean(social.dat$social.diff.scale[social.dat$treatment.assignment == 0])
outmat['social.perceptions','int.sd'] = sd(social.dat$social.diff.scale[social.dat$treatment.assignment == 0])
outmat['social.perceptions','seg.mean'] = mean(social.dat$social.diff.scale[social.dat$treatment.assignment == 1])
outmat['social.perceptions','seg.sd'] = sd(social.dat$social.diff.scale[social.dat$treatment.assignment == 1])
outmat['social.perceptions','ate'] = ate
outmat['social.perceptions','se'] = sd(distout.ci)
outmat['social.perceptions','ci.low'] = quantile(distout.ci,.025)
outmat['social.perceptions','ci.high'] = quantile(distout.ci,.975)
outmat['social.perceptions','p.value'] = disp$greater.p.value
outmat['social.perceptions','cohens.d'] = cohen$estimate
outmat['social.perceptions','N'] = nrow(social.dat)



###physical analysis
diff.Xs = as.matrix(physical.dat[,x.vars])
diff.ri.perms = genperms(physical.dat$treatment.assignment,
                         blockvar = physical.dat$Block.ID)
diff.ri.probs = genprob(diff.ri.perms)
ate = estate(Y = physical.dat$physical.diff.scale,
             Z= physical.dat$treatment.assignment,
             prob = diff.ri.probs,
             X = diff.Xs)
diff.Ys = genouts(Y = physical.dat$physical.diff.scale,
                  Z= physical.dat$treatment.assignment,
                  ate = 0)
distout <- gendist(diff.Ys,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
diff.Ys.ci = genouts(Y = physical.dat$physical.diff.scale,
                  Z= physical.dat$treatment.assignment,
                  ate = ate)
distout.ci <- gendist(diff.Ys.ci,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = physical.dat$physical.diff.scale,
                f = as.factor(physical.dat$treatment.assignment))


outmat['physical.perceptions','int.mean'] = mean(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 0])
outmat['physical.perceptions','int.sd'] = sd(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 0])
outmat['physical.perceptions','seg.mean'] = mean(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 1])
outmat['physical.perceptions','seg.sd'] = sd(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 1])
outmat['physical.perceptions','ate'] = ate
outmat['physical.perceptions','se'] = sd(distout.ci)
outmat['physical.perceptions','ci.low'] = quantile(distout.ci,.025)
outmat['physical.perceptions','ci.high'] = quantile(distout.ci,.975)
outmat['physical.perceptions','p.value'] = disp$greater.p.value
outmat['physical.perceptions','cohens.d'] = cohen$estimate
outmat['physical.perceptions','N'] = nrow(physical.dat)


##format the matrix for output
outmat = as.data.frame(round(outmat,2))
outmat[,1] = paste(outmat[,1], ' (',outmat[,2],')',sep='')
outmat[,3] = paste(outmat[,3], ' (',outmat[,4],')',sep='')
outmat[,5] = paste(outmat[,5], ' (',outmat[,6],')',sep='')
outmat[,7] = paste('[',outmat[,7],',',outmat[,8],']',sep='')

outmat[,8] = NULL; outmat[,6] = NULL; outmat[,4] = NULL; outmat[,2] = NULL

outmat[,'cohens.d'] = outmat[,'cohens.d']*(-1)


##output Appendix Table 2
outmat.reduced = outmat[,c('ate','ci.low','p.value','cohens.d','N')]
table.2 = xtable(outmat.reduced)
rownames(table.2) = c('Allocation','Attribution', "Social Perceptions", "Physical Perceptions")
colnames(table.2) = c('Beta (SE)','CI','P','D','N')

print(table.2,
      "Table2.tex",
      type = 'latex',
      floating = F,
      booktabs = T)


##output Appendix Table S3
table.s3 = xtable(outmat,
                  digits = 3)
rownames(table.s3) = c('Allocation','Attribution', "Social Perceptions", "Physical Perceptions")
colnames(table.s3) = c('Integrated (SD)','Segregated (SD)','Beta (SE)','CI','P','D','N')

print(table.s3,'TableS3.tex',
      type = 'latex',
      floating = F,
      booktabs = T)


cat('\n \n Experiment 2 Results: \n')
print(outmat.reduced)


##Table S4, Balance Table
out.balance.test =  xBalance(fmla = treatment.assignment ~  male + non.white + hispanic + college.grad + conservative  + age + income.recode  + weight + self.height,
                             data = all.dat, 
                             report = c("std.diffs","z.scores","chisquare.test"), 
                             strata = factor(all.dat$Block.ID),
                             na.rm=T)

balance.test.results <- out.balance.test$results[,c("std.diff","z"),]   

use.vars = c('treatment.assignment','male', 'non.white', 'hispanic', 'college.grad', 'conservative', 'age', 'income.recode', 'weight','self.height')
mean.dat = all.dat[,use.vars]
mean.dat = mean.dat[complete.cases(mean.dat),]
mean.table = t(aggregate(mean.dat,by = list(mean.dat$treatment.assignment),
                         FUN = mean,na.rm=T,
                         simplify = T)) 
mean.table = mean.table[3:nrow(mean.table),]
bind.table = cbind(mean.table,balance.test.results)
bind.table = rbind(bind.table,c(nrow(all.dat[all.dat$treatment.assignment == 0,]),nrow(all.dat[all.dat$treatment.assignment == 1,]),NA,NA))
row.names(bind.table) = c("Male","Non White", "Hispanic", "College Graduate","Conservative", "Age","Income","Weight", "Height","N")
colnames(bind.table) = c("Integrated","Segregated","Std Diff", "Z-score")
table.s4 = xtable(bind.table,
                   digits = c(2,2,2,3,3))
print(table.s4,'TableS4.tex',
      type = 'latex',
      floating = F,
      booktabs = T)
